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Abstract. Ultracold polar molecules, in highly anisotropic traps and interacting via a repulsive 
dipolar potential, may form onc-dimensional chains at high densities. According to classical 
theory, at low temperatures there exists a critical value of the density at which a second order 
phase transition from a linear to a zigzag chain occurs. We study the effect of thermal and 
quantum fluctuations on these self-organized structures using classical and quantum Monte Carlo 
methods, by means of which we evaluate the pair correlation function and the static structure 
factor. Depending on the parameters, these functions exhibit properties typical of a crystalline 
or of a liquid system. We compare the thermal and the quantum results, identifying analogies 
and differences. Finally, we discuss experimental parameter regimes where the effects of quantum 
fluctuations on the linear - zigzag transition can be observed. 
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1. Introduction 

The experimental realization of strongly interacting atomic ensembles at ultralow 
temperatures allows the observation of self-organization into ordered patterns. 
Prominent examples are Coulomb crystals of laser-cooled ions in Paul and Penning 
traps [lj. In these systems laser cooling permits one to access the crystalline regime, 
which emerges from the competition of the long-range Coulomb interaction and the 
trapping potential at low temperature. Different structures are found depending on the 
trap geometry and the number of atoms [H El El S] - In one of the early experimental 
works [2], [3] , the phase diagram of the various low dimensional crystalline structures was 
determined as a function of the trapping frequency of the harmonic confinement and of 
the number of atoms, starting from a setup with the ions arranged in a linear chain. 
Later works, concerned about the possibility to measure quantum fluctuations affecting 
the structure of these systems, pointed out that the ground state of ion crystals is 
essentially classical [5] . Recent theoretical works, nevertheless, predicted that for chains 
of few ions quantum fluctuations, such as single-particle tunneling at the structural 
transition linear - zigzag, may be observed [6] . 

Under this perspective, ultracold dipolar gases of atoms or polar molecules are 
interesting systems, as the nature of the dipolar interaction allows one to observe 
the interplay between quantum degeneracy and long-range forces. Dipolar gases of 
ultracold atoms have been experimentally realized with Chromium atoms [8]. In 
these experiments the strength of the s-wave scattering interaction is conveniently 
controlled by tuning the Feshbach resonance and can be made vanishing, thus leading 
to the realization of purely dipolar systems. Stability against collapse is warranted by 
polarizing the gas in a two-dimensional geometry, such that the dipolar interactions 
can be considered as purely repulsive 0, [10]. Several theoretical works have studied 
the phases of dipolar gases of atoms or polar molecules as a function of density and 
dimensionality [Til [121 KH1 H31 US1 LTS] - In the context of this paper, particularly relevant 
are studies which predicted the creation of self-organized structures in a two-dimensional 
bulk system, showing that the ground state may exhibit the typical features of a crystal 
or of a quantum fluid depending on the density (9[ [TO] . 

One-dimensional dipolar gases can be experimentally realized in highly anisotropic 
traps [IT]. In this limit, Luttinger liquid models describe the long-range properties of 
dipolar gases at ultralow temperatures [H] . In this case, by tuning the density, the phase 
of the gas undergoes a crossover from a Tonks-Girardeau gas to a crystal-like phase [19|. 
Differing from Coulomb systems, the crystal-like phase is energetically favorable at large 
densities, as it is typical of short-range interacting systems. 

In Ref. [20] we presented a theoretical study of low dimensional dipolar gases 
focusing on the effects of quantum fluctuations on the formation of transverse 
correlations at the classical critical point from a linear to a planar distribution, where 
the one-dimensional structure is unstable with respect to increasing the atomic density 
and/or to decreasing the transverse potential. In the present article, we provide further 
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details of that study, and extend it to considering thermal effects in the formation of 
quasi-crystalline structures in the classical domain. The analysis of the thermal and 
quantum fluctuations is performed numerically, by means of Monte Carlo simulations. 
The pair correlation function and the structure form factor for the linear and zigzag 
chains for the classical case at finite temperature and the quantum case at T = are 
compared, and analogies and differences discussed. 

This article is organized as follows. In Sec. [2] we introduce the theoretical model for 
a two-dimensional gas of dipolar bosons in presence of tight transverse confinement. In 
Sec. [3] we discuss the mechanical stability of a one-dimensional array of dipoles, assuming 
classical dynamics and neglecting thermal fluctuations. In Sec. H] the numerical methods 
applied for simulating the thermal and quantum fluctuations are described. The effect 
of thermal fluctuations on the quasi- crystalline phase are discussed in Sec. [5j The phase 
diagram for the quantum ground state is introduced and discussed in Sec. [6] while the 
numerical results for the pair correlation function and structure form factor of linear 
and zigzag phases are presented in Sec. [71 The conclusions are drawn in Sec. El and 
details of the calculations are reported in the appendices. 



2. The model 



We consider N bosonic atoms or molecules of mass m at ultralow temperatures. The 
particles have induced or permanent dipole and are confined by an external harmonic 
potential on the plane x — y at z = by means of a very steep confinement along the 
z axis, while their dipoles are aligned along the z direction by means of an external 
field, such that their mutual interaction is always repulsive. By appropriately tuning 
an external magnetic field to the top of a Feshbach resonance, the s-wave scattering 
can be made negligible so that the particles interact exclusively through the repulsive 
long-range dipolar interaction potential, which here reads 

Vint = , | CM | 3 , (1) 

4vr|p! - p 2 | 3 

where Cdd is the dipolar interaction strength and pi = (a^,^) is the particle position 
vector in the x — y plane. 

The Hamiltonian of the dipoles contains the kinetic energy, the harmonic 
confinement and the dipole-dipole interaction 



H = y ZL + l mu 2 y z +±yy Ldd - (2) 

3 1 J 3=1 j¥* lH HJI 

where Pj is the conjugate momentum of the coordinate Pj and v t is the confinement 
frequency of the particles along the ^-direction. In the following we will consider a 
homogeneous system along the x-direction with periodic boundary conditions. For 
the classical regime, where the particles are distinguishable, this will correspond to 
imposing pi + N = Pi- In the quantum regime, this consists in imposing that the 
ground-state wave function fulfills the condition ip(xi, . . . , Xj + L, . . . , xn, Vi ■ ■ ■ , Vn) = 
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ip(xi, . . . ,Xj, . . . , xat, yi, ... , y N ) for all j, where L is the size of the system along the 
x-direction. 



3. Mechanical stability of a classical chain of dipoles 

In this section the dipolar system is studied in the classical regime. We first consider the 
equilibrium positions of the dipoles as determined by the trap and the repulsive dipolar 
potential, and compute the critical value of the density at which a harmonic expansion 
around these equilibrium positions is unstable. 

For later convenience, we rescale the Hamiltonian in terms of a unit energy for the 
transverse oscillator, Eho, and define the rescaled classical variables 7ij = pj / \J mEh 



and Qj = pj/ \J Ehojmv'l \ where now Qj 
Hamiltonian H = H/ Eho reads 



(xj,yj). In these units, the dimensionless 



(3) 



where the potential energy V is given by 



j 



E 



and r = r /y E ho /mi>? with 



mu?C, 



r 



t^dd 



(4) 



(5) 



a characteristic length, whose physical meaning will become more transparent when 
considering the quantum mechanical case. The classical equilibrium positions of 
potential are here denoted by the vector gf = (xf\y\ ), which is solution of 
the equations 



dV 

dxi 
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We assume the interparticle distance along the x axis to be uniform, and given by the 
length 

a = l/n, (8) 



where n is the linear density and we denote by n = Eho/mu^ the rescaled density. 
Within the assumption of equispaced dipoles, Eq. §6§ is automatically satisfied. In order 



to find a solution of Eq. ([7]), we assume the ansatz gy = (ja, (— 1)- ? 6/2), with a = 1/n, 
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Figure 1. (color online) Transverse variance of the dipolar distribution for the classical case, 
\J (y 2 ), (in units of a/, = \J Eho/rnv 2 ) as a function of the linear density fi. The variance 
is evaluated from the classical equilibrium positions of the potential, and it does not contain 
the thermal fluctuation: it hence gives the displacement of the dipoles from the trap axis at 
equilibrium. The value of b in classical system can be directly compared to the mean transverse 
spreading y/ (y 2 ) as both are expected to coincide in linear and zigzag chains at low temperatures. 
The critical value n c where the transition linear-zigzag chain occurs is shown by the arrow. The 
solid line corresponds to the solution of Eq. ([9]) , the symbols to the results of a classical Monte 
Carlo simulation, the dashed line is a Taylor expansion of Eq. ([9]), see Eq. (|A.10|) . Here, nf = 3. 



and which for 6 = corresponds to a linear arrangement of the dipoles, while for 6^0 
it gives a zigzag structure. Using this ansatz, Eq. ([7]) becomes 



b 
2 



N/2 



12f 



{ [(2£-l) 2 5 2 + 6 2 ] 5 /2 



0. 



(9) 



which relates the value of b to a. Notice that b = is always solution of Eq. Q. Values 
of b 7^ are solutions of Eq. (Q provided that the linear density h is larger than some 



critical value n > n c where 



£r, 



=-1/5 



o (10) 

is the critical density which depends on the trap frequency v t , £ = (8/(93 ^(5))Y^ 5 = 
0.6078 ... is a constant, and ( is the Riemann zeta function. The values of b, computed 
by solving numerically Eq. (Q for different values of n, are shown in Fig. [1] and compared 
to values of a/ (y 2 ) obtained by classical Monte Carlo simulations (see Sec. H] for details 
on the numerical simulations). 

The dispersion relations of the dipolar chain inside a harmonic potential have been 
determined in Ref. [22]. In this work we consider a simpler situation, assuming there 
is no harmonic potential along the x direction, and periodic boundary conditions are 
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imposed at the edges. In this regime the frequencies of the phononic modes of the 
harmonic crystal are simply found using a phonon-wave ansatz, see for instance [21] [23] . 
This ansatz solves exactly the equations of motion for the longitudinal and transverse 
displacements from the equilibrium positions, q~j = Xj — ja and Wj = jjj, respectively, 
which read 

= - 12f °Eife-^ ( n ) 



1 



fij = -*i + 3f: o^-^(^-*i+ i ) (12) 

m 3 

One finds two branches, for the axial and transverse excitations, which are given by 

2 fu \ 48f 1 • 2 k 3 a f-\ q\ 

II fl ^1 „o 2 

j>0 J 

w±(k) = l- — ^-sm— (14) 
i>o J 

where k = 2irl/Na with I = 0, 1, 2, . . . , N — 1. The two branches are analogous to an 
acoustic (cuy) and an optical (uj±) branch. A few relevant modes are depicted in Fig. [2j 
The acoustic branch, Eq. ( Tl~3|) . is characterized by a sound velocity c = 
A/l2f C(3)/a 3 . This value is consistent with the value found from the relation [24] 

mc 2 — ndfi/dn, (15) 

where \x is the chemical potential, \i = dE/dN and E is the ground state energy, which 
at large densities is equal to E = (A^ 4 / L 3 )((3)Cdd/4n (see |Appendix A| and Ref. [19]). 



From the dispersion relation of the optical branch, Eq. (Tl4l) . one finds the stability 
condition of the linear chain configuration. In fact, Eq. ([14]) gives imaginary values of 
the transverse frequency for values of the linear density n < n c , with h = 1/a and n c 
the critical density. Its value is found by solving the relation minfc(u;_i_) = 0, which gives 
n c in Eq. fflUj) and k — ir/a. The corresponding mode whose frequency becomes zero at 
h c is the zigzag mode (see Fig. [2b), 

^soft = (_ 1)^/2, (16) 

where b is the amplitude of the mode when the chain is stable. When one can assume 
that the dipolar chain exhibits long range order (at very large interaction strengths), the 
structural instability at n = n c can be characterized as a classical second order phase 
transition, applying a straightforward extension to dipolar interactions of the analysis 
presented in Ref. [21] for Coulomb interactions. The details are reported in Appendix 
A. 



4. Monte Carlo methods 



The numerical study of the properties of a many-body system is carried out using 
Monte Carlo methods. The advantage of these methods is that they permit to evaluate 
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Figure 2. Some eigenmodes of the linear chain and of the zigzag structure, a) Linear chain: 
long- wavelength axial mode at wave vector k = 2ir/Na, b) short-wavelength transverse mode at 
wave vector k = ir/a — 2tt/Ncl, and (c) transverse mode at k = ir/a (zigzag mode), d) Zigzag 
structure: long wavelength eigenmode (this mode can be obtained as a superposition of the modes 
of the linear chain shown in a) and b)). 



multidimensional integrals over the phase space. In a classical system of dimensionality 
D this correspond to an integral over 2 x N x D variables R = {ri, r^; vi, . . . , vjv} 
with Tj and Vj being position and velocity of a particle j, with j = 1, . . . , N. In a 
quantum system the dimensionality of the phase space is iV x D with R = {ri, r^} 
in coordinate representation. The average of some quantity of interest A(R) has the 
form 

j ... j A(R)p(K) dR 



(A) 



(17) 



J...J P (R) dR 

with p(R) the averaging weight. In a classical system at thermal equilibrium position 
and velocity are independent variables. For a quantity A(R) = A(ti, ...,rjv) depending 
only on particle coordinates, the velocities in Eq. (ITTj) can be integrated out. As a result 
the average of interest is performed by only integrating over the position variables, using 
the averaging weight 

V(r 1 ,...,r N ) 



cxp 



18) 
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where V is the potential and T the temperature. 

In the quantum case the averages are performed using the coordinate representation. 
In a variational (VMC) calculation the probability distribution is defined by the modulus 
squared of the trial wave function ip T {ri, ...,y n ), 

p V Mc(ri,-..,r N ) = |?/; T (ri, ...,rAr)| 2 . (19) 

The trial wave function ip? depends usually on some variational parameters, which 
can be optimized by minimizing the variational energy of the system. Once the 
optimal parameters are known, this wave function is used as guiding wave function 
in Diffusion Monte Carlo (DMC) calculations. The DMC method solves the imaginary- 
time Schrodinger equation for the product of the wave function ip and the guiding wave 
function ipT- For large values of imaginary times the DMC method permits one to 
sample the probability distribution 

PDMc(ri, -, rjv) = ^r(ri, rjv)0 o (ri, r N ), (20) 

where 0o( r i> •••> t n) is the ground state wave function (more details on the DMC method 
can be found, for example, in Ref. [25J). 

In both cases, classical (1181) and quantum ( 1191) . the probability distribution is known 
analytically and can be written explicitly. This permits one to perform its sampling 
using the Metropolis algorithm. A trial move converting the probability from p to p' 
is always accepted if p' > p (in a classical system this means that energy in the new 
configuration has a lower energy or, being the same, it has a larger value of the wave 
function in the quantum case), otherwise the move is accepted with probability p'/p. We 
generate a new configuration by displacing one particle by a random distance, which we 
draw from a Gaussian distribution. The width of the Gaussian distribution is adjusted 
so that on average we have 50% acceptance rate. Under those conditions Metropolis 
algorithm guarantees that the generated set of configurations (Markov chain) has the 
desired distribution and the average (A) is stochastically evaluated as an average over 
the Markov chain. 

We use the Bijl-Jastrow construction of the trial wave function 

N N 

m*i, -> ^) = n n Mi* - r *D- ( 2i ) 

i=l j<k 

The one-body term is chosen in Gaussian form /i(r) = exp{— ay 2 /a 2 ho \. In the limit of 
quasi-one-dimensional system, where the excitations of the radial harmonic confinement 
levels are absent, the optimal value a = 1/2 corresponds to the ground state wave 
function of a harmonic oscillator. If the system leaves the quasi-one-dimensional regime, 
the system spreads in the radial direction, and accordingly the value of a is reduced. 
The two-body term is the same as used in Refs. [TQl [20] 

( dKopy/Daho/r), < r < i? par 

f 2 (r) = 1 C 2 exp{-C 3 a ho [l/r + l/(L x - r)]}, R par <r<L/2 (22) 
[ 1, L/2<r , 
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with K (r) being the modified Bessel function of the second kind. Coefficients C±, C 2 , C 3 
are fixed by the conditions of the continuity of the wave function and its first derivative. 
The parameter < R par < L/2 is free and is fixed by a variational procedure (the 
dependence of the energy on this parameter is rather weak and in many cases the 
choice Rp ar = L/4 is sufficient). When the distance between two particles is small, the 
influence of other particles can be neglected. As a result, /2( r ) is well approximated by 
the solution of the two-body scattering problem (see short distance behavior of (122]) ). 
This makes the calculation to be stable even for a divergent interaction potential. We 
symmetrize the argument of the long-range part as 1/r + 1/(L X — r). In this way we 
ensure the zero derivative of the trial wave function at the boundary r = L x /2 according 
to the periodic boundary conditions. 

DMC averaging ( l20j) provides exact estimators for operators commuting with the 
Hamiltonian. In particular, the ground state energy is found in a statistically exact 
way. Instead, density profiles and other functions of particle positions are obtained as 
"mixed" estimators (A)dmc = ( , 0t|^-|^o)/('0t|0o) which are biased by a particular 
choice of the trial wave function. Ideally, one has to find a pure estimator, i.e. 
averages over the ground state (A) = (0 o |^4|0o)/(0o|0o)- The optimized trial wave 
function is expected to be close to the true ground state wave function, so that the 
difference \5ip) = \i/jt) — \4>o) is small and one can use the mixed estimator together 
with the variational one (A)vmc = (V'tI^I^t) / (V'tIV't) (obtained according to (1T9"]) ) to 
extrapolate the pure estimator. The two extrapolation procedures 



have the same order of accuracy and contains errors which are quadratic in \5ijj)- 
Alternatively, for local quantities one can use more involved techniques of pure 
estimators [26J. For simplicity we will use the rules ( 12311241) and check that both 
estimators agree within the accuracy of interest. 

5. Thermal fluctuations around the classical equilibrium positions 

The structures discussed in Sec. |3] were determined by simply evaluating the classical 
equilibrium positions of the potential energy, Eq. (jll), as a function of the linear density 
and of the transverse trapping frequency. In this case, we identified a sharp condition, 

_ 1/5 a k ~ ~ 

n c = Sr , as given by Eq. (1101) . such that for n < n c the equilibrium positions 
are along a line (linear chain), otherwise are arranged along multiple lines. For an 
interval n c < h < Um, with um a fixed value, the structure is a zigzag chain. Thermal 
fluctuations will smear the sharp transition at n c . In this section we compute their 
effect around n c for various values of the temperature using a classical Monte Carlo 
calculation. 

For later convenience, we report the temperature in units of the energy Eh a , which 
we set Eho = hu t . In order to study the transition linear - zigzag structure we evaluate 



(A) 
(A) 



2(A)daic — (A)vmc 
( a )dmc/ (A)vmc 



(23) 
(24) 




Figure 3. (color online) Contour plot of the mean transverse size of the system (y 2 ) at 
T = lO~ e hv t /kB as a function of the parameters fo and n. The contour lines correspond to 
fixed values of (y 2 ) (in units of a\ ), the scale is reported on the right side of the plot. The dashed 
line, obtained from solving Eq. (|10j) . separates the linear from the zigzag structures; the solid line 
separates the zigzag from the multiple chain structures. The calculations are made for a system 
of N — 20 particles, where we covered the plane (h - fo) with a grid of 50 x 50 size, each point 
corresponding to a run with a different choice of parameters. Effects of the finite grid are visible 
in Monte Carlo data in the vicinity of the dashed line. 



the transverse variance (y 2 ) and report the numerical results as a function of the linear 
density n and the interaction strength fo in a contour plot for two different values of the 
temperature, see Figs. (3JH1 In addition, in the figures we plot two lines corresponding 
to the transition from a linear to a zigzag (dashed line), and from a zigzag to a 
multiple chain (solid line), as obtained by finding the classical equilibrium positions 
of the potential (j4j), and plotted for comparison. 

Figure [3] reports the result obtained setting T = \Q~ & hv t /kB- Here, one can still 
identify a sharp transition from linear to zigzag, and from zigzag to multiple chains, 
according to the predictions given by Eq. [101 Such situations could be found in setups 
with extremely high transverse frequencies, although it is questionable whether such 
trap could be experimentally realized. Figure H] displays the case T = Tu/f/fes, showing 
that at this temperature the transition is completely smeared out. 

We now study how thermal effects modify the pair correlation function, which is 
related to the probability of finding two particles separated by a distance (x,y), and is 
here defined as 

g 2 (x,y) = (Af( Pl )Af(p2)), (25) 

where M(p) = Ylf &{p ~ Pi) is the two-dimensional density, and g2(x,y) depends only 
on the relative distances x = x\ — x 2 and y — yi — yi- Here, integrating out the radial 




coordinate from the two-dimensional density gives the linear density j M{p) dy = n(x) 
and its average is equal to the linear density (n(x)) = n. 

The pair correlation function is displayed in Fig. [5] for various values of the 
temperature. Here, the "hole" for x = and y = is due to the hard-core repulsive 
interaction potential, which forbids two particles being at the same position. There is a 
clear strong correlation between neighboring particles, which is lost as the temperature 
and the interparticle distance increase. This is a manifestation of the absence of a true 
crystalline order. 

Information over the properties of the system can be inferred by means of Bragg 
spectroscopy, giving access to the static structure factor. The latter is given by 

S(k) = ±Wf- 1l ) (26) 

where A/k = / e tkp N{p)dp = ,=i exp{zkpj} is the Fourier transform of the two- 
dimensional density N(p). The right column of plots in Fig. [5] displays the static 
structure factor, corresponding to the pair correlation function plotted on the left side. 
Here, one sees the appearance of a high peak, at k x = 2ir/a, and of a second peak, at 
k x = Air/a, at low temperatures. The peaks would correspond to Bragg scattering by 
an ordered structure at interparticle spacing a, and are thus signatures of very strong 
correlations in the system. Their height diminishes, while they become broader, as 
temperature increases and the axial structure is lost. The disappearance of the second 
peak at higher temperatures is typical of a gas phase. 

Figure [6] displays the pair correlation function and the structure form factor in the 
parameter regime, in which the zigzag structure is expected at ultralow temperatures. 
As before, the plots of g 2 (x,y) show a hole for the zero separation, caused by the 
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Figure 5. (color online) Pair distribution function (n(f?)n(0)} (left) and static structure factor 
(Sfc) for the linear structures (right). Temperatures are (rows from upper to low) T = 0.05, 
T = 0.1, T = 0.2 and T = 0.5hv t /k B . Number of particles is N = 20, f = 7.5, n = 0.4 



dipolar repulsion. Moreover, one finds a crystal-like ordering at very low temperatures, 
T m 0.05hu t /kB, which disappears as the temperature is increased by a factor of 10. 
In addition, at low temperatures one observes large probability to find two particles 
with the same transversal displacement, such that yi — yj = 0, with doubled periodicity 
with respect to the linear chain. One also observes additional peaks corresponding to 
correlations in the y direction, corresponding to the doubled chain (zigzag) structure. 
As in one half of the cases the particle of reference belongs to the upper (lower) chain, 
there are peaks in g2(x,y) for y = —2b (y = 2b). The peaks at y = have larger 
amplitudes as they are sum of the contributions from both upper and lower chains. 
The static structure factor is plotted as a function of the two-dimensional wave 
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Figure 6. (color online) Same as in Fig. but with fo = 100 and n = 0.4, warranting the 
zigzag regime. Temperatures are (rows from upper to low) T = 0.05, T = 0.1, T = 0.2 and 
T = O.bhvt/kB- Differing from Fig. [5l the static form factor is now a mesh plot as a function of 
h x (in units of 2n/a) and of k y (in units of 2ir/b). 

vector (k x , k y ). We note that different boundary conditions in x and y lead to different 
accessible values of moment umj|| At low temperatures the static structure factor reminds 
the geometry of a two-dimensional crystal, and naturally shows how a one-dimensional 
structure transforms into a two-dimensional one. There is a number of directions in 
which well-defined peaks are observed. The peaks in direction (k x ,k y ) = (ir/a, 0) are 

| Here, periodic boundary conditions along x lead to a discretization in momentum, with unit step 
equal to 2tt/L = (2i:n)/N , while in y direction the momentum is continuum. In practice this means that 
typically the width of the crystal-like peaks along k x cannot be determined, while along k y the peaks 
have well-defined width. In the thermodynamic limit the discretization disappears and the momentum 
is continuous in all directions. 
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Figure 7. (color online) Static structure factor of a linear chain of polar molecules. Parameters 
are the same is in Fig.[SJ while the temperatures are T = 10~ 4 ; 3x 10~ 4 ; 10~ 3 ; 3x 10~ 3 ; 10~ 2 hv t /kB 
(see plot). 

similar to the ones shown in Fig. [5] for the linear chain. The peaks in the transverse 
direction (k x ,k y ) = (0, ir/b) or in diagonal directions (k x ,k y ) = (n/a,ii/b), etc. are 
signatures of the realization of a zigzag structure. As the temperature is increased, the 
number of peaks is decreased until only the central peak is left. We note, however, that 
the presence of the central peak is somehow trivial as it counts the number of particles 
in the system. 

Finally, we study how the ordering in the linear chain is destroyed as the 
temperature is increased from T = 10~ 3 hut/ks, up to T = 0.1hu t /kB, by looking at 
the features of the structure form factor. In Fig. [7] is plotted as a function of k x 
for N = 20 particles, for various values of the temperature. The peak at k = counts 
the number of particles and its height equals to 20 in all cases. The height of the other 
peaks, at multiples of 2n/a, diminishes as T is increased. In particular, it decreases with 
the number of the peak. Moreover, the height of the peaks seems to follow a power-law 
dependence as a function of the peak number. These observations will be compared 
with the effect of quantum fluctuations on ordering at T = 0. 

6. The quantum ground state 

The definition of a one-dimensional ground state density distribution of dipoles in the 
quantum regime requires the existence of an energy gap for exciting the transverse 
motion, such that at T = the relevant excitations of the dipolar gas are essentially 
along the longitudinal axis. In order to quantify this statement we first introduce the 
characteristic lengths of the system. The confinement in the transverse direction sets 
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the characteristic length (oscillator length) 

a ho = \fh I mv t , (27) 

which determines the width of the single-particle ground-state wave packet, and which 
we choose as unit length. We extend the classical model of Eq. ([3]) by taking the unit 
energy Eh = Kv t in which now the classical variables q and 7r are replaced by the 
corresponding quantum conjugate operators. The potential energy V is given by the 
Eq. (j3J). For a one-dimensional configuration, when the transverse excitations are frozen 
out, the parameter r gives the characteristic length scale of interactions, and has to 
be compared to the length scale of quantum fluctuations set by n -1 , where n is the 
linear density. The system is in a gas or quasi-crystal phase depending on whether the 
product nr$ is much smaller or much larger than unity. For nr$ ^> 1 the dipolar gas 
is in the quasi-ordered phase with ground-state energy E£ D) = N(nr ) 3 ((3)h 2 /(mr 2 ), 
which is the potential energy of a classical crystal [19] (see also Appendix A ). In the 
quantum gas regime, for nr <C 1, the properties of the system are well described by 
the Tonks-Girardeau model [27] with ground state energy = N7T 2 h 2 n 2 /(Qm) (the 

expression for e!£q does not depend on the strength of the dipolar interaction, as in 
this regime the interparticle distance is large and the potential interaction enters only 
through the positions of the nodes of the ground state wave function, while the potential 
energy can be neglected) . These two regimes are shown in the phase diagram in Fig. [HI 
displaying the various phases of the ground state of the dipolar system as a function of 
f and n = naho, and are separated by the grey short-dashed line, which indicates the 
curve nro = nro = 1- 

Dynamically the system is one-dimensional when the chemical potential /i is much 
smaller than the level spacing of the transverse oscillator, /i <C hu t . In this regime, the 
quantum state of the system can be described within the Luttinger liquid formalism, as 
it was shown in Ref. [18], and no true long-range order is found at finite densities [28J. 
The curve E^ 1D ^/N = fiv t separates the one- and two-dimensional phases and is shown 
in Fig. [8] by the black dashed line. In the quasi-ordered phase, for nr ^> 1, this 



corresponds to the inequality E&' /N <C Hv t , and which leads to the relation n<Cr 
In the quantum gas regime, for nr^ <C 1, the condition of being one-dimensiona 
E^q^ /N <C Kvt, which is equivalent to the requirement naho *C 1. 



The phase diagram in Fig.[8]can be experimentally explored by changing the density 
n, moving vertically, and by varying the transverse frequency u t , moving parallel to the 
straight short-dashed line. Typical parameters of experiments with one-dimensional 
gases are given for instance in Ref. [T7J: for a^o ~ 35 nm, 60 — 400 atoms per tube of 
length varying between 15 and 50/xm, one finds f ~ 0.07, h pa 0.04 — 0.3. Among all 
atom species with which condensation has been reached, 52 Cr has the largest value of 
r = 2.4 nm [5J. Larger values of r can be reached using polar molecules such as CO, 
ND3, HCN, CsCl with ro = 5 nm — 340 /xm, permitting one to cover regions of the 
phase diagram up to the classical region (nr 3> 1). 



Thermal and quantum fluctuations in chains of ultracold polar molecules 



16 



n 



10" 



10- 1 



2D 



Tqnks-Girardeau 
quasi 1D 



i i i 1 1 1 1 1 — i — i i ■ ■ i ■ | 1 — i — i i ■ ■ ■ ■ | r 




10"' 



io- 1 



10" 



10 1 



10' 



10 J 



10 4 



Figure 8. (color online) Ground-state phase diagram as a function of the parameters fo and 
n. The dashed black line corresponds to the curve E 1D /N = hv t and identifies the two regimes, 
where the dynamics is essentially one- or two-dimensional. The Tonks-Girardeau and the classical 
crystal limits are explicitly indicated in the plot (dashed-dotted lines). The short-dashed grey line, 
at nr = 1, separates gas- and solid-like phases. These lines are to be intended as indicators, the 
transition of the gas from one phase to another being a crossover. The green solid lines separate 
different quasi-crystalline structures in the classical system. Inset: portion of the phase diagram in 
the solid-like regime and close to the region, where the classical transition linear-zigzag structure is 
found (solid line) . The red circles are results of quantum Monte Carlo calculations and correspond 
to the appearance of a double-peak structure in the radial density profile (see details in Ref. [20]). 
The size of the symbols denotes the error bars. 



7. Effects of quantum fluctuations on ordering 

In this section we present numerical results for the pair correlation function and the 
structure form factor of the quantum gas of bosons at T = 0, focusing on the solid-like 
regime and studying in particular the parameter regimes, where classically the linear 
and zigzag chains are observed. Now, the structural properties are defined from the 
interplay between the quantum kinetic term, dipolar interactions and the confinement 
potential in the Hamiltonian ([3]). The results are compared with the classical case at 
finite temperature, discussed in Sec. [5j 

Figure [9] displays the pair correlation function and the static structure factor of the 
quantum system at zero temperature. For the same parameter regimes in a classical 
system, we observed at very low temperature quasi- crystalline linear structures, see 
Fig. [5] The results in Fig. [9] hence clearly show that quantum fluctuations significantly 
modify the phase of the system. No features of quasi-ordering can be observed: the 
pair correlation function is almost uniform (apart for the hole at x = y = 0), and the 

f As discussed in Sec. the hole comes from the short-range divergence of the repulsive potential. This 
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g 2 (x, y) 




Figure 9. (color online) Pair distribution function g 2 (x,y) (left) and static structure form factor 
Sk for the quantum gas of dipolar bosons at T = 0, with N = 20 particles, n = 0.4 and fo = 7.5. 
The results are obtained with a Quantum Monte Carlo simulation, see Sec. SJ Extrapolation (f?5|> 
is used to remove bias due to a particular choice of trial wave function 



static structure form factor exhibits a small peak at the wavevector, corresponding to 
the mean linear density n, while no further peaks at multiples of k = 2im are visible. An 
analogous situation is observed in the classical system at sufficiently high temperatures, 
compare with Fig. [51 

Figure [10] displays the pair correlation function and the static structure factor of 
the same quantum system at zero temperature for the parameter regime, in which one 
observes a quasi-ordered zigzag structure in the classical system, compare with Fig. El 
Now, quasi-ordering of the particle is visible in the pair correlation function, which 
corresponds to the appearance of well-defined peaks in the static structure form factor. 
Nevertheless, one does not observe secondary peaks, showing that there is still no true 
crystalline structure. Comparison with the classical case, Fig. [6J shows similarities to 
the finite-temperature behavior in a regime when there is no ordering, while strong 
correlations are still present. 

Let us now comment on the behavior of the static structure form factor in the 
quantum regime. In a dipolar quantum system the low-momentum part is described 
by phonons. The static structure factor grows linearly with k for small momenta (see 
Fig. ED: 

^) = i ^ ^ 

where c is the speed of sound, see also Eq. (fT5l) . At zero temperature in the k — > 
limit the static structure factor vanishes S(k) — > C0- The low-momentum behavior is 
important as it reflects the number fluctuations in a large system ((5N) 2 ) = NS(0) as it 

can be easily understood by considering the lower limit of the integral giving the potential energy of 
dipolar interactions: / V(r)g2(r) dr oc J <?2(r)/l r l 3 dr. This energy diverges if 32(f) does not vanish at 
short distances. Consequently the pair correlation function has to vanish when two particles meet both 
in classical and quantum system. In a quantum system this also means that the two-body Jastrow term 
has to vanish at short distances, which is in agreement with the behavior of the two-body scattering 
solution used in the construction of the trial wave function (|22p . 

I Note that from the definition (|26|) the value of the S(k) is discontinuous at zero momentum where it 
gives the number of particles for k — 0. For a gas S(k) — NS(k) is instead continuous. 
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Figure 10. (color online) Same as in Fig. [9J but with fo = 100. The static structure form factor 
is displayed as in Fig. [6] 



is discussed in details in Appendix B refer to Eq. (IB. 31) . On the other hand, the number 
fluctuations at finite temperature are related to the isothermal compressibility [29] 
4 = -V 2 /(mN) (dP/dV) T as 

This means that at finite temperatures S(k) no longer vanishes in the zero momentum 
limit, instead it goes to a constant value S(k) — > ksT/mc^, k — > 0, defined by the 
temperature and the compressibility ( fT5l) . Moreover, the effects of the temperature in 
quantum Bose systems on the low-momentum part of the static structure factor can be 
found analytically [30] and S(k) is given by the following expression 

~/,s Kk , hkc knT h 2 k 2 .,. , m ;x _ . . 

S « = 2^ Ct W = ^ + IW + -' < 30 > 
The low-momenta dependence is no longer linear, but a quadratic one with the coefficient 
of proportionality entirely defined by the temperature. 

The one- dimensional case has been studied in the literature [HI [28j EI], and 
evidence has been given that the one-dimensional dipolar systems manifest Luttinger 
liquid behavior even if the interactions are of a long-range type [18]. From Haldane's 
approach [32] it is possible to write asymptotic expressions for the one-dimensional pair 
correlation function g2(x). Although this effective description is valid only at large 
distances, it is sufficient to predict the behavior of the peaks in the static structure 
factor using properties of the Fourier transformation (1261) . The static structure factor 
has been computed for various values of hr in Refs. [28| [3Tj . showing that at large 
values multiple peaks appear. Their height decays with the peak order according to a 
power law decay (with the exponent depending on the Luttinger parameter, and hence 
on the linear density [31]), showing interesting analogies with the classical case at finite 
temperatures. Clear important differences between quantum mechanical and thermal 
effects in the structure form factor can be identified in its behavior at low values of k, 
as the value of the static structure factor at zero momentum is related to the integral of 
the density fluctuations over the whole volume of the system. In the quantum mechanical 
case, it grows from zero linearly in k (see Fig. [9]) while in the classical case at finite T 
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no longer vanishes at zero momentum, instead it goes to a constant value, defined by 
the temperature and compressibility (see Appendix B for details). This can be lead 
back to the issue how number fluctuations {(SN) 2 ) increase with the volume size L D , D 
being dimensionality of the homogeneous space. Indeed, thermal fluctuations increase 
linearly with the volume, while fluctuations in quantum systems at zero temperature 
follow the law ((5N) 2 ) oc L D ~ l InCL [331 El], where C is a constant. In particular, in a 
one-dimensional system (and also in a quasi-one- dimensional geometry), the fluctuations 
follow the logarithmic dependence {(5N) 2 ) oc InL [35] . 

8. Conclusions 

We have discussed the properties of the classical and quantum distribution of a low- 
dimensional gas of dipolar particles at thermal equilibrium, by computing numerically 
the pair correlation function and the static structure form factor, focusing to the regime 
where one would expect quasi-ordered structures of dipoles in a linear or zigzag structure. 
Thermal and quantum fluctuations tend to destroy ordering, bringing further evidence 
that a true self-organized crystal cannot exist at finite temperature or (in a quantum 
system) at a finite interaction strength. Analogies and differences between the effects of 
thermal and quantum fluctuations on the phase diagram, the pair correlation function 
and the static structure factor are discussed. Behavior of the static structure factor 
for small momenta is addressed. In the limit of large dipolar interaction strengths or 
of large linear density, the transition from a one-dimensional to a planar configuration 
occurs with the creation of a mesoscopic structure in the transverse direction, which 
exhibits the main features of the transition from a single to a double and then a multiple 
chain distribution of dipoles, while quasi-order is observed also in the transverse pair 
correlation function at large densities. Such patterns are characterized by non-local 
correlations, which arise from zero-point fluctuations, and which may be important 
resources for the realization of quantum simulators [361 E7] . 

Support by the ESF (EUROQUAM "CMMC"), the European Commission 
(EMALI, MRTN-CT-2006-035369; SCALA, Contract No. 015714) and the Spanish 
Ministerio de Educacion y Ciencia (Consolider Ingenio 2010 "QOIT"; FIS2007-66944; 
FIS2008-04403; Ramon-y-Cajal; Juan de la Cierva) are acknowledged. 

Appendix A. Transition from linear chain to the zigzag structure 

For very large interaction strengths, the linear chain becomes unstable for n > n c . 
Under this condition one can show that, within a certain interval of densities, the stable 
configuration is a planar zigzag structure. The zigzag structure is composed by two 
chains, one shifted with respect to the other by a lattice distance a, so that with respect 
to the linear chain the periodicity doubles from a to 2a. Following [2Tj, we here make 
the Landau theory of the classical phase transition from a linear chain of dipoles to a 
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zigzag, which is strictly valid when thermal effects can be neglected. 

We identify the order parameter of the zigzag as the distance 6/2 of the dipoles 
from the trap axis x and the control parameter as the density n while the transverse trap 
frequency v t remains fixed. A similar analysis can be performed using the transverse 
frequency as the control parameter and fixing the density, see [2T] . In order to study the 
behavior close to the instability point, we expand Eq. (TJJ up to the fourth order around 
the equilibrium positions of the chain, V = Ylt=o where I labels the order. Here, 

is the ground state energy in the limit of strong correlations, which for N ^> 1 
reads 

y(°)=JV^. (A.l) 
a 6 

Moreover, as we are expanding close to the equilibrium positions of the linear chain, 
the first order term V^ l > vanishes. Using the decomposition into the eigenmodes of the 
linear chain, the quadratic term takes the form 

^ (2) = ^EE(-nW 20 i s)2 +^)^ )2 ) ( A - 2 ) 

k s=± 

where 

OA 1 12f ° 1 ' 2 kia /a q \ 

i>0 



and it coincides with ujj_(k) 2 , Eq. (JED, for n < h c , while the normal modes 0^ and 

(s) 

ty k are related to axial and transverse displacements by the relation 

& = V^J2{ e k +) cos kja + Q^ sin fcja) (A.4) 



W 3 



k 

^ E cos k i a + sin k i a ) ( A - 5 ) 

k 



with k = 2-Kl/Na with I — 0,1,2, ... , N/2 and the subscript ± indicates the parity by 
the reflection k — > —k. For convenience, we denote by \l/o the mode of the linear chain at 
wave vector ko = ir/a (zigzag mode). The linear chain becomes mechanically unstable 
when, by varying h, the frequency of the zigzag mode with wave vector k = ir/a, 
Eq. (Tl6l) . vanishes. Around the instability point the zigzag mode is significantly coupled 
to other quasi-degenerate modes by the third and fourth order terms and 
These quasi-degenerate modes are long wavelengths axial modes Q$k at wave vectors 
5k, such that \5k\a <C 1, and short wavelength transverse modes ^ko+Sk at wave vector 
k = ko + 5k, with \5k\a <C 1. At lowest order in the expansion in 5ka, the effective 
potential for the modes which are relevantly coupled at the instability reads 

Sk>0 

where for brevity we denote ask = a{ko — 5k) and = ^k +8k- The fourth order 
potential for the relevant modes is given by: 

5k>0 
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where 



^=^fl-^K(7) (A. 



N V 2? , 

and the function / contains a sum of products of three amplitudes tygk f° r 8k 0) 
whose precise form is not important for the present discussion. We now allow the linear 
density n to take values in the interval [n c — 5, n c + 5] , such that ask may take only 
small but negative values, and for a>sk < we determine the corrections to the 
equilibrium positions of the linear chain using Eq. ( 1A.6j) . assuming that the role of the 
fourth order potential V k ~ ko is to give rise to a small displacement b <C a with respect 
to the equilibrium interparticle distance a. Using a stability analysis of the zigzag 
configuration in a similar way as the one used for the linear-zigzag transition in low 
dimensional ion crystals [21], we find that the stable solution has ^sk = for 5k > 
and 

where «o = 1 — (n/n c ) 5 < 0. Correspondingly, the spatial displacement b reads 
b = 2^q/\/N and, using Eqs. ffTU]) and Eq. flA.81) it can be written as 



b = Ba\\l-^ (A.10) 
where we have introduced the dimensionless constant 



B= .l*m = 0.634... (A.11) 
Y 5 127C(7) V ; 

Note that the expression for b found in Eq. (1A.10I) corresponds to the second order 
expansion of Eq. (Q in the parameter b/a. The exact result found by solving numerically 
Eq. ([9]) and the approximated result Eq. ( 1A.10|) are compared in Fig. [TJ 



The jump in the classical ground-state energy at the transition point is given by 
the difference between the potential of the soft mode at each sides of the critical points. 
We denote by 

V soit = ^l + A^ (A. 12) 



the effective potential for the soft mode, where A is given by Eq. ( 1A.8I) . At zero 
temperature, for n = n c — 5, with 5 > then V soit = 0, while for n = n c + 5 then 
V soit ^ 0. Letting 5 -> 0+ we find that the energy difference per particle is given by 

AE = ^ - ni)-V"*(n n;) = _ , )2 ^ 

where C = [4650C(5)]/[127C(7)] = 37.65. . .. Notice that the second derivative of AE 
with respect to n is discontinuous at the critical point. At zero temperature this means 
that the phase transition is of the second order. 



Thermal and quantum fluctuations in chains of ultracold polar molecules 



22 



Appendix B. Static structure form factor and thermodynamic properties 

Small momentum behavior is directly related to the density fluctuations in the system. 
Indeed, the number fluctuations in volume L is related to the integral of pair correlation 
function 

{(8N) 2 ) = / dx\ I dx 2 (5n(xi) 5n(x 2 )) = N + / dxi / dx 2 [g 2 (xi — x 2 ) — n 2 } (B.l) 

JL JL JL JL 

Here we consider the case of a linear chain considering it as a one-dimensional objects 
so that the volume is linear L and the pair correlation function saturates to a constant 
value, g 2 {x) — > n 2 as \x\ — > oo. The static structure factor is related by the Fourier 
transformation to g 2 (x) according to relation g 2 {x) = n 2 + n j(S(k) — l)e~* fc:!: dk/2n, 
which is inverse to Eq. (1261) . By substituting this expression to Eq. ( IB. II) one gets 

((5N) 2 ) = N + n j dx 1 J^dx 2 J ° ^(S(k) - l)e~ lkx (B.2) 

Changing to center of mass variables x = X\ — x 2 ; X = (x\ + x 2 )/2 one separates 
spatial integrals into two. One is trivial J L dX = L. The second gives the ^-function 
J e lkx dx = 2n5(k) leading to a simple relation 

((5N) 2 ) = NS(0). (B.3) 
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